#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>

#ifndef _WIN32
#include <sys/time.h>
#include <sys/sysinfo.h>
#include <unistd.h>
#else
#include <windows.h>
#endif

typedef struct {
        double xL, fL, xR, fR, area0, tol;
} package_t;

int max_number_of_threads;
/* Global storage for the thead count and a protective mutex: */
int thread_count = 0;
pthread_mutex_t thread_count_mutex = PTHREAD_MUTEX_INITIALIZER;

/* Global storage for the integral and a protective mutex: */
double sum = 0.0;
pthread_mutex_t sum_mutex = PTHREAD_MUTEX_INITIALIZER;
int thread_id = 0;

/* The function to integrate: */
double f(double x){
        return exp(-x*x);
}

void * rec_quad(void *argIn){

	/* Local variables: */
	double tol, xL, xR, fL, fR, area0;
	double xM, fM, aL, aR;
	pthread_t *threadL, *threadR;
	package_t *lPackage, *rPackage;

	package_t *inPck = (package_t *)argIn;
	/* Extract into local variables: */
	tol = inPck->tol;
	xL = inPck->xL; xR = inPck->xR;
	fL = inPck->fL; fR = inPck->fR;
	area0 = inPck->area0;

	xM = 0.5*(xL + xR); fM = f(xM);
	aL = 0.5*(xM - xL)*(fM + fL); aR = 0.5*(xR - xM)*(fM + fR);

	if (fabs(area0 - (aL + aR)) / area0 <= tol)
	{
		pthread_mutex_lock(&sum_mutex);
		sum += aL + aR;
		pthread_mutex_unlock(&sum_mutex);
	}
	else
	{
		pthread_mutex_lock(&thread_count_mutex);
		if (thread_count + 2 <= max_number_of_threads)
		{ /* Do work with threads... */
			thread_count += 2;
			pthread_mutex_unlock(&thread_count_mutex);
			/* Create two threads to do this subintegral: */
			lPackage = (package_t *)malloc(sizeof(package_t));
			if (!lPackage)
			{
				fprintf(stderr, "No memory\n");
				exit(3);
			}
			else
			{
				lPackage->tol = tol;
				lPackage->xL = xL; lPackage->fL = fL;
				lPackage->xR = xM; lPackage->fR = fM;
				lPackage->area0 = aL;
			}
			rPackage = (package_t *)malloc(sizeof(package_t));
			if (!rPackage)
			{
				fprintf(stderr, "No memory\n");
				exit(3);
			}
			else
			{
				rPackage->tol = tol;
				rPackage->xL = xM; rPackage->fL = fM;
				rPackage->xR = xR; rPackage->fR = fR;
				rPackage->area0 = aR;
			}
			threadL = (pthread_t *)malloc(sizeof(pthread_t));
			/* Spawn two new threads to do each subinterval: */
			if (pthread_create(threadL, NULL, rec_quad, (void*)lPackage) != 0)
			{
				perror("pthread_create");
				exit(2);
			}
			threadR = (pthread_t *)malloc(sizeof(pthread_t));
			if (pthread_create(threadR, NULL, rec_quad, (void*)rPackage) != 0)
			{
				perror("pthread_create");
				exit(2);
			}
			/* Wait for the threads to join and clear used memory */
			pthread_join(*threadL, NULL);
			free(lPackage); free(threadL);
			pthread_join(*threadR, NULL);
			free(rPackage); free(threadR);
		}

		else
		{
			pthread_mutex_unlock(&thread_count_mutex);

			double threadSum, dx, xL_i, xR_i;
			int niters, sc;

			threadSum = 0.0;
			niters = (int)(((xR - xL) / sqrt(tol)) + 1);
			dx = (xR - xL) / niters;
			for (sc = 0; sc < niters; sc++){
				xL_i = xL + dx*sc;
				xR_i = xL_i + dx;
				threadSum += 0.5*dx*(f(xL_i) + f(xR_i));
			}
			pthread_mutex_lock(&sum_mutex);
			sum += threadSum;
			pthread_mutex_unlock(&sum_mutex);
		}
	}
}

void test_with_time(void *istruct){
		
	double _duration;
#ifndef WIN32
	struct timeval tim1,tim2;   
    gettimeofday(&tim1, NULL);   
#else
	LARGE_INTEGER frequency;
	LARGE_INTEGER t1,t2;
	QueryPerformanceFrequency(&frequency);
	QueryPerformanceCounter(&t1);
#endif
        rec_quad(istruct);		
#ifndef WIN32	
	gettimeofday(&tim2, NULL);
	_duration=tim2.tv_sec+(tim2.tv_usec/1000000.0)-tim1.tv_sec+(tim1.tv_usec/1000000.0); 
#else
	QueryPerformanceCounter(&t2);
	_duration=(double)(t2.QuadPart-t1.QuadPart)/frequency.QuadPart;
#endif

	printf("Integral = %18.10f\n", sum);
	printf("Time = %18.10f\n",_duration);
}

int main(){
	package_t *istruct;

	double xL = -4.0;
	double xR = 4.0;
	double tol = 0.003;

	double fL = f(xL);
	double fR = f(xR);
	double area0 = 0.5*(xR - xL)*(fL + fR);	

	max_number_of_threads = 20;

	/* Package this information into a package type: */
	istruct = (package_t *)malloc(sizeof(package_t));
	if (!istruct)
	{
		fprintf(stderr, "No memory\n");
		exit(2);
    }
	
	istruct->tol = tol;
	istruct->xL = xL; istruct->fL = fL;
	istruct->xR = xR; istruct->fR = fR;
	istruct->area0 = area0;
	/*1*/
	test_with_time(istruct);
	/*2*/
	tol = 0.00003;
	istruct->tol = tol;
	sum = 0;
	test_with_time(istruct);
	/*3*/
	tol = 0.0000003;
	istruct->tol = tol;
	sum = 0;
	test_with_time(istruct);
}